1.1. Surface Matching : example on human scapula


This notebook presents a method to match a point cloud taken with any scanning system to a digital 3D scan in stl format.

1.1.1. Introduction

1.1.1.1. Method high lines

Surface matching is based on an optimization strategy. The objective of this optimization is to identify the parameters of the reference change allowing to pass from the coordinate system of the STL file to the coordinate system of the scan system.

This change of reference frame is formalized by a passage matrix formed by a rotation vector and a translation vector according to the Rodrigues formalism. There are thus a total of 6 DOF.

The optimization parameters of the cost function will therefore be these two vectors. The residual is defined by the signed distance between the points coming from the measurement systems after the change of reference frame and the stl mesh. This distance is calculated between the points of the closest scanning system and the STL mesh.

This problem being badly formulated, it is necessary to establish an initial solution in order to have a correct result.

This is why in a first step, we establish an initial using points that are identified on the STL model but also on the real model. A minimum of 3 points are necessary to block the 6DOF.

Then, in a second step, an optimization will be performed with the whole scan taking as initial solution the one previously calculated.

1.1.1.2. Modules importations

import trimesh
import numpy as np
from scipy import optimize
import plot_mesh as pltm
import plotly.graph_objects as go
import pandas as pd
import gbu

1.1.1.3. Case of study

Left Human scapula stl file

fig = pltm.create_mesh3D(stl_file="scapula.stl", title="Scapula Left <br> stl file")
fig.show()

1.1.2. Define reference points on STL file

We define reference points on the STl that can be easily found on the real model A point on the glenoid, coracoid and acromion is selected.

p_glene_stl = np.array([0.02515289854664963, 0.03122873596800735, 0.036798809060995745])
p_coracoide_stl = np.array(
    [0.026787610816186087, 0.05609095153197788, -0.0043762069034760384]
)
p_acromion_stl = np.array(
    [-0.008420164259823573, 0.07140986464971243, 0.018059532550748165]
)
fig = pltm.add_points(points=p_glene_stl, fig=fig, name="glène")
fig = pltm.add_points(points=p_coracoide_stl, fig=fig, name="coracoide")
fig = pltm.add_points(points=p_acromion_stl, fig=fig, name="acromion")
fig.show()

The same points are easily found on the model

../_images/scapula.jpg

1.1.3. Computations

1.1.3.1. Get scanned points from our system scan device

df_scan = pd.read_pickle("data_scan.p")
df_scan
p_glene_scapula p_acromion_scapula p_coracoide_scapula p_scan_scapula
p3d p3d p3d p3d
x y z x y z x y z x y z
0 0.014046 0.010801 0.118422 -0.023197 -0.022039 0.086035 -0.002085 -0.037214 0.123242 0.015003 -0.004692 0.11943
1 0.014152 0.01073 0.118276 -0.023097 -0.022086 0.085978 -0.002091 -0.037263 0.123189 0.015018 -0.004673 0.119397
2 0.014185 0.010854 0.11838 -0.023238 -0.022144 0.086024 -0.002113 -0.037188 0.123366 0.015072 -0.004562 0.119343
3 0.014362 0.01094 0.118554 -0.023217 -0.022149 0.086317 -0.001946 -0.037163 0.123109 0.015048 -0.004463 0.119347
4 0.014379 0.010899 0.118494 -0.023224 -0.022023 0.086214 -0.001959 -0.036706 0.12319 0.015013 -0.004507 0.11929
... ... ... ... ... ... ... ... ... ... ... ... ...
79 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.015566 -0.001602 0.126646
80 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.01562 -0.001514 0.126536
81 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.015541 -0.001624 0.126567
82 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.015546 -0.001497 0.126532
83 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.015578 -0.001497 0.126524

84 rows × 12 columns

1.1.3.2. Plot scanned points in their reference frame

Here the points are expressed in scapula reference frame

cols = set([cols[0] for cols in df_scan.columns])
fig = go.Figure()
for i, c in enumerate(cols):
    fig.add_trace(
        go.Scatter3d(
            x=df_scan[c].p3d.x,
            y=df_scan[c].p3d.y,
            z=df_scan[c].p3d.z,
            name=c,
            marker_size=2,
            mode="markers",
        )
    )
fig.show()

1.1.4. Create initial guess

1.1.4.1. Gathering input data of reference points

p_glene_scan = np.array(df_scan["p_glene_scapula"].p3d.mean())
p_coracoide_scan = np.array(df_scan["p_coracoide_scapula"].p3d.mean())
p_acromion_scan = np.array(df_scan["p_acromion_scapula"].p3d.mean())

tri_point_stl = np.concatenate([p_glene_stl, p_coracoide_stl, p_acromion_stl]).reshape(
    -1, 3
)
tri_point_scan = np.concatenate(
    [p_glene_scan, p_coracoide_scan, p_acromion_scan]
).reshape(-1, 3)

1.1.4.2. Compute initial guess

def unpack(X):
    return X.reshape(2, 3)


def cost(X):
    Rsc2stl, Tsc2stl = unpack(X)
    tri_point_out = gbu.utils.apply_RT(tri_point_scan, Rsc2stl, Tsc2stl)
    dist = abs(tri_point_stl - tri_point_out)
    pen_pin_radius = 0
    res = dist - pen_pin_radius

    return res.flatten()


X0 = np.zeros(6)
sol = optimize.least_squares(
    cost,
    X0,
    method="lm",
    ftol=1.0e-12,
    xtol=1.0e-12,
    gtol=1.0e-10,
)

X_pre = sol.x
R_inital_guess, T_initial_guess = unpack(X_pre)
print(f"Initial guess solution: rvec={R_inital_guess}, tvec={T_initial_guess}")
Initial guess solution: rvec=[ 2.23440972 -0.14576349  0.55625112], tvec=[-0.02428643  0.12368749  0.09225085]

1.2. Surface match with full scan

1.2.1. Gathering input data of full scan

p_full_scan_sc = np.array(df_scan["p_scan_scapula"])

# load full scapula expand
scale = 1e-3  # convert to meter
mesh_scapula = trimesh.exchange.load.load("scapula.stl", type="stl")
mesh_scapula.vertices *= scale

1.2.2. Optimization

def cost(X, p3d, mesh):
    Rsc2stl, Tsc2stl = unpack(X)
    p_stl = gbu.utils.apply_RT(p3d, Rsc2stl, Tsc2stl)
    dist = trimesh.proximity.signed_distance(mesh, p_stl)
    pen_pin_radius = 0
    res = dist - pen_pin_radius
    return res


sol = optimize.least_squares(
    cost,
    X_pre,
    method="lm",
    ftol=1.0e-12,
    xtol=1.0e-12,
    gtol=1.0e-10,
    args=(p_full_scan_sc, mesh_scapula),
)

R_sc2stl_sol, T_sc2stl_sol = unpack(sol.x)
print(f"Final solution: rvec={R_sc2stl_sol}, tvec={T_sc2stl_sol}")
Final solution: rvec=[ 2.2283109  -0.16742782  0.55011184], tvec=[-0.02341616  0.12482507  0.09038642]

1.2.3. Plot results

p_full_scan_stl = gbu.utils.apply_RT(p_full_scan_sc, R_sc2stl_sol, T_sc2stl_sol)

fig = pltm.create_mesh3D(stl_file="scapula.stl", title="Scapula Left <br> stl file")

fig = pltm.add_points(
    points=p_full_scan_stl, fig=fig, name="scanned_points", marker_size=10
)

fig.show()